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We study the collapse of rapidly rotating supermassive stars that may have formed in the early 
Universe. By self-consistent ly simulating the dynamics from the onset of collapse using three- 
dimensional general-relativistic hydrodynamics with fully dynamical spacetime evolution, we show 
that seed perturbations in the progenitor can lead to the formation of a system of two high-spin 
supermassive black holes, which inspiral and merge under the emission of powerful gravitational 
radiation that could be observed at redshifts z > 10 with the DECIGO or Big Bang Observer 
gravitational- wave observatories, assuming supermassive stars in the mass range 10 4 — 10 6 M©. The 
remnant is rapidly spinning with dimensionless spin a* = 0.9. The surrounding accretion disk 
contains ~ 10% of the initial mass. 

PACS numbers: 04.25.D-, 04.30.Db, 97.60.Bw, 02.70.Bf 



The observation of supermassive black holes (SMBHs) 
with masses M > 1O 9 M in the early Universe at red- 
shifts z > 7 (e.g., [1 ) calls for new pathways for explain- 
ing the existence of such massive black holes when the 
Universe was less than 1 Gyr old (see, e.g., [2] [3] for re- 
cent reviews). Since the time available for the formation 
of a massive black hole via accretion from an initial seed 
black hole is quite short, it appears likely that the seed 
objects were quite massive themselves. Sufficient initial 
seed mass could be offered by the theoretical possibility 
of supermassive stars (SMSs) with masses 10 4 — 1O 6 M . 
Since SMSs can become general-relativistically unstable 
[U |5], they can collapse and either form supermassive 
seed black holes, or, perhaps, die in powerful thermonu- 
clear explosions with energies > 10 55 ergs [6HH]- 

A popular scenario for the formation of an SMS is the 
direct collapse of a primordial gas cloud in dark matter 
halos with virial temperatures > 10 4 K [9HT6]. The re- 
quired rapid accretion rates [17] can be maintained by 
avoiding early fragmentation, either by isothermal col- 
lapse in the absence of H2 cooling [lOj [T8H22] . or by 
turbulent accretion [T4HT6] . Pulsational instabilities or 
radiative feedback, which could limit the accretion rate, 
are sufficiently small to pose no problem for rapid growth 
[23j [24] . Since the primordial gas cloud is likely to carry 
substantial angular momentum, the centrifugal barrier 
must be overcome during the collapse of the halo gas, 
which could be achieved by angular momentum transfer 
via gravitational torques [T4] . 

In [25], the collapse of a uniformly rotating SMS was 
studied in axisymmetry. Later, [26j [27] investigated the 
three-dimensional collapse of differentially rotating SMSs 
to SMBHs and found that the collapse proceeds axisym- 



metrically. These studies were complemented by [28J [29] , 
who investigated the collapse of rapidly differentially ro- 
tating SMSs with small seed perturbations that lead to 
fragmentation during collapse. These authors found that 
a small initial m = 1 density perturbation grows expo- 
nentially and leads to a single, off-centered fragment that 
collapses to a black hole. In the case of an initial m = 2 
density perturbation, two orbiting and collapsing frag- 
ments form. Before the two fragments form black holes, 
however, a run-away collapse at the center of the SMS 
occurs, leading to a single black hole at the center of the 
star. These studies used a T = 4/3 T-law equation of 
state (EOS), which is appropriate for SMSs dominated 
by radiation pressure. 

Recently, [6] investigated the collapse of SMSs in ax- 
isymmetry using a microphysical EOS that models the ef- 
fects of radiation and electron-positron pair production. 
They found that above central temperatures > 10 9 K, 
pair creation effectively reduces the adiabatic index be- 
low T = 4/3, thus accelerating the collapse once these 
temperatures are reached. 

In this Letter, we show that by assuming a T-law EOS 
with slightly reduced T < 4/3, e.g., because of pair- 
creation, it is possible to form a binary black hole system, 
provided certain small seed perturbations are present at 
the onset of collapse and the star is rapidly differentially 
rotating. The binary black hole system subsequently in- 
spirals and merges under the emission of powerful gravi- 
tational radiation. 

Methods. We employ general-relativistic hydrodynam- 
ics with adaptive mesh-refinement provided by the open 
source EinsteinToolkit [30] [31] using the Llama mul- 
tipatch infrastructure [32] [33] , and a modified version of 



TABLE I: Black hole Christodoulou masses Mbh and dimen- 
sionless Kerr spin parameter a BH for each black hole in each 
model. In model M2G2, we list the parameters of the inspi- 
raling two black holes in the first two rows, and the merger 
remnant by the end of the simulation in the third row. We 
also report the disk mass Mdisk by the end of the simula- 
tion, the measured accretion rate M, and the radiated GW 
energy Eqw- For model M2G2, where multiple resolutions 
are available, we are able to provide error bars. Units are in 
c— G — M — K — 1, unless otherwise specified. 



M1G1 



M2G1 M2G2 



BH mass Mbh [M] 
BH spin a BH 



5.4 



0.9 



5.8 



0.9 



3.0 ±0.1 
3.0 ±0.1 
5.8 ±0.2 
0.7 ±0.02 
0.7 ±0.02 
0.9 ±0.01 



bar. disk mass Mdisk [M] 1.3 1 


0.7 ±0.2 


accretion rate M 1.2 x 10~ 3 2 x 10~' 


4 6.7 x 10 -5 


rad. GW energy E GW [%] 0.02 0.16 


3.71 



WEN05 reconstruction [30] 153] . We evolve the space- 
time geometry with the Z4c system [35], extract gravi- 
tational waves (GWs) at future null infinity via Cauchy- 
characteristic extraction [36H39] , and find apparent hori- 
zons using AHFinderDirect [40]. We employ a T-law 
EOS and use an artificial low-density atmosphere (10 -10 
the central density). If not stated otherwise, all units 
are m. c = G = M = K = 1, where K is the polytropic 
scale. Since our models can be rescaled to any desired 
mass, we report most numbers in units of mass M. For 
convenience, we give the conversion factors to cgs units 
in the Appendix. 

Initial conditions. We consider 3 initial SMS mod- 
els. Their initial data are generated via Hachisu's self- 
consistent field method [41] [42], which requires as in- 
put the central density p c of the star, and a polar-to- 
equatorial axes ratio r p /r e between and 1 together with 
a dimensionless parameter A to set the degree of differen- 
tial rotation. All models are set up as marginally stable 
T = 4/3 polytropes with scale K = 1. In every case, 
the central density is set to p c = 3.38 x 10 -6 M -2 , the 
axes ratio is set to r p /r e = 0.24, and the differential 
rotation parameter is set to A = 1/3. With these in- 
put parameters, the baryonic mass of the star becomes 
M*,bar = 7.0527M, the gravitational mass of the star 
is M* 5 adm = 7.0037M, and the angular momentum be- 
comes J = 52.206 M -2 , which corresponds to a dimen- 
sionless Kerr spin parameter of a* = J/M 2 = 1.0643. 
The resulting equatorial radius is r e = 82M. Due to the 
rapid rotation, the star has quasi-toroidal shape (its den- 
sity maximum is off-center) . The various models differ by 
the type of initial perturbation that is applied and by how 
the collapse is induced. In two models, M1G1 and M2G1, 
the collapse is accelerated by reducing the polytropic 
scale K by 10 -3 initially. During evolution, the two mod- 



els use an adiabatic index Y = 4/3. In model M2G2, the 
collapse is triggered by reducing T = 4/3 to T = 1.33. To 
induce fragmentation during collapse, we apply a density 
perturbation to the initial configuration. The perturba- 
tion is given by pi n i — »> pi n i (1 ± A m r sin(m</>)), where 
TTi > is an integer, A m is the perturbation ampli- 
tude, r a cylindrical coordinate radius, and <\> a cylin- 
drical coordinate angle. This perturbation offers reason- 
able overlap with the corresponding quasinormal modes 
of the star. In model M1G1, we apply a m = 1 per- 
turbation, and in models M2G1 and M2G2, we apply a 
m = 2 perturbation. We use a perturbation amplitude 
of A m = 10- 3 /r e « 1.22 x 1(T 5 . 

Grid setup. We use the multiblock grid setup de- 
scribed in [32]: A central Cartesian grid is surrounded 
by six spherical "inflated-cube" grids. Their interface is 
located at radius R s = 240M, and the outer boundary 
is at Rb = 4000M. The central Cartesian grid is capa- 
ble of 2:1 mesh refinement. We use 2 additional finer 
grids, one covering the entire star, the other covering the 
high density torus. We perform simulations using three 
resolutions labeled by rO, rl, and r2. In our baseline 
resolution rl, the initially finest level has a resolution 
of Ax = 0.4M. During evolution, we track each form- 
ing fragment with a moving refinement center. While a 
fragment collapses, we progressively add an additional 
refinement level every time the fragment's central den- 
sity increases by another order of magnitude. In total, 
we switch on up to 3 additional refinement levels for each 
fragment until a black hole forms. The finest level has 
a resolution of Ax = 0.05M. The outer spherical grids 
have radial spacing Ar = 3.2M and use N ang = 21 cells 
per patch and angular direction, which corresponds to 
an angular resolution of ~ 4.3°. The other resolutions 
rO and r2 have 25% decreased and 25% increased resolu- 
tion. See the Appendix for demonstration of numerical 
convergence. 

Dynamics. Following the initial pressure reduction, 
either by reducing K or by reducing T, the rapidly ro- 
tating SMSs undergo accelerated collapse. Depending on 
the initial density perturbation, the stars fragment, and 
collapse accelerates in the fragments. Model M1G1, due 
to its m = 1 initial density perturbation, forms a single 
off-centered fragment which orbits around the star's cen- 
ter while collapsing. At T = 2145 M, an apparent horizon 
emerges around the center of the fragment. Model M2G1, 
due to its m = 2 initial density perturbation, forms two 
orbiting fragments that slowly inspiral and collapse. Af- 
ter T ~ 2335M, when the two fragments are still well 
separated, run-away collapse occurs at the center of the 
star, leading to a rapid increase of the central density. 
A single apparent horizon emerges around the center of 
the SMS at T = 2470M. Models M1G1 and M2G1 were 
both considered in [28, 29 , and we confirm their find- 
ings, though we are able to continue the evolution well 
beyond the appearance of an apparent horizon. We stop 
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FIG. 1: Snapshots of a slice of the equatorial density distri- 
bution of model M2G2. Dark colors indicate high density, 
light colors indicate low density. The logarithmic density col- 
ormap ranges from 10" 7 M~ 2 (white) to 10" 3 M~ 2 (black). In 
the bottom two panels, the colormap is rescaled to the range 
[10~ 8 M~ 2 , 10~ 4 M~ 2 ] for the sake of presentation. The up- 
per two and the bottom right panels show an extent of ±40M, 
while the remaining panels show an extent of ±20M. 



the simulation at T ~ 3500M. Model M2G2 initially fol- 
lows qualitatively the evolution of model M2G1. Due to 
its softer T, however, the two fragments collapse much 
faster, forming a pair of two black holes. In Fig. [I] we 
show the equatorial density distribution of model M2G2 
at various stages during the evolution. The upper left 
panel shows the inner 80M of the initial density distri- 
bution. The toroidal high density ring is clearly visible. 
The upper right panel shows a snapshot during collapse. 
Due to the initial density perturbation, a strong m = 2 
deformation arises with two inspiraling high-density frag- 
ments. The center-left panel shows the density distri- 
bution at time T = 1130M when an apparent horizon 
appears within each of the two still well separated frag- 
ments. The apparent horizons are indicated by white 
circular regions. For the purpose of visualization, they 
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FIG. 2: Top panel: Evolution of the maximum density for all 
models. Center panel: Spin and mass evolution of all three 
black hole horizons of model M2G2. Lower panel: (£,m) = 
(2, 2) spherical harmonic mode of the "+" polarization of the 
emitted gravitational radiation rescaled by distance D. 



are approximated by the lapse gauge function a (with 
the region a < 0.19 closely resembling their coordinate 
shapes). The two nascent black holes have identical mass 
and spin. Initially, they have a coordinate separation 
D = 13.2M , a Christodoulou mass M BH = 0.899M, and 
dimensionless spin ag H = 0.286. Their spin axes are 
both aligned with the orbital angular momentum. Trans- 
fer of angular momentum from the differentially rotating 
torus onto the black holes via accretion leads to an out- 
wards spiraling motion until the two black holes reach 
a separation D = 24.7 M. From there, they start to 
inspiral driven by GW emission and merge after ~ 1.5 
orbits. The tracks of the black holes are indicated by 
red (solid and dashed) lines in the center-left panel of 
Fig. [I] In the center-right panel, the black holes have 
completed close to one orbit. Around each black hole, 
material is dragged and accreted, forming spiral patterns. 
During inspiral, while further accreting material, each of 
the two black holes grows to a mass Mbh = 3 ± 0.1M 
and dimensionless spin parameter <2g H = 0.7 ± 0.02 just 
before merger (center panel of Fig. [2]). In the bottom 
left panel of Fig. [T] the two black holes are about to 
merge at time T = 1680M, and a common apparent 
horizon appears. The bottom-right panel shows the sit- 



uation at the end of the simulation after the system has 
settled to a quasi-stationary state. The merger rem- 
nant has a Christodoulou mass Mbh = 5.8 ± 0.2M, 
and dimensionless spin ag H = 0.9 ± 0.01. The sur- 
rounding accretion disk quickly settles to a baryonic mass 
M disk = 0.7 ± 0.2M. 

GW emission. All models emit GWs during collapse 
and while the fragments form and orbit around the cen- 
ter of the star. The dominant radiated GW mode for 
all models is (£,m) = (2,2). In models with m = 2 
deformation, the emission is particularly strong. The 
strongest emission is generated by model M2G2 due to 
the two inspiraling and merging SMBHs. In Fig. |2j we 
show the evolution of the maximum density as well as 
the (£, m) = (2,2) mode of the "+" polarization of the 
emitted GW signal. As the maximum density increases, 
the oscillatory GW signal rises in amplitude. Follow- 
ing BH formation, a GW ring-down signal is emitted in 
models M1G1 and M2G1. In model M2G2, following 
BH formation, we obtain a binary black hole inspiral sig- 
nal which increases in frequency and amplitude towards 
merger. Following merger, the remnant emits ring-down 
GWs. The entire GW signal is in the dHz frequency 
band for M* ~ 10 4 M SMSs, and in the mHz band for 
M* - 10 6 M SMSs. 

Using modes up to £ = 8, we compute the radiated 
energy emitted in GWs. In model M1G1, 0.02% of the 
initial gravitational mass M^adm is radiated in GWs, 
0.16% is radiated in model M2G1, and 3.71% is radi- 
ated in model M2G2 (see Table [[]). We compute the 
detectabilities by the proposed space-borne eLISA, DE- 
CIGO, and Big Bang Observer GW observatories [43]- 
[45] over a range of redshifts z € [5, 100] and source 
masses M* e [1O 4 M , 1O 6 M ]. Assuming a ACDM cos- 
mology with parameters measured by the Planck mission 
[46] , we compute the luminosity distance D(z) and red- 
shifted mass (1 + z)M* for a given redshift z. The op- 
timal and angle- averaged signal-to-noise ratios (SNRs) 
(e.g. [47] [48]) can then be obtained from the Fourier- 
transformed angle-dependent GW signal h(6,(j)) using 
the corresponding theoretical GW observatory sensitivity 
curves [43] [45j [49] , the luminosity D(z) and correspond- 
ing redshift ed mass (1 + z)M*. Assuming a minimum 
SNR of 8 for detection, we obtain a maximum redshift of 
z ~ 25 for a M* = 10 4 M SMS in DECIGO and the Big 
Bang Observer, provided the system's spin axis is point- 
ing towards Earth ((0,0) = (0,0)). Both detectors are 
limited by the white dwarf confusion noise at low frequen- 
cies (e.g. [49 ) and have very similar sensitivity for the 
considered mass range. Higher mass SMSs are detectable 
only out to smaller redshifts. For a M* = 10 6 M SMS, 
we obtain a maximum redshift z c± 16 for detection. By 
assuming a random orientation of the spin axis, we ob- 
tain a mean GW detectability out to z c± 23 for low-mass 
SMSs, and z ~ 13 for high-mass SMSs. In eLISA, the 
signal is barely detectable at z ~ 6, which is outside the 



relevant range of redshifts z > 10 at which SMSs are 
anticipated to exist. 

Discussion. We have self-consistently simulated the 
collapse of SMSs using 3 + 1 general-relativistic hydro- 
dynamics with approximation-free dynamical spacetime 
evolution. We have shown that it is possible to form 
a coalescing binary SMBH system in SMS collapse pro- 
vided the following conditions are realized: (i) the SMS 
must be rapidly differentially rotating, (ii) the gas pres- 
sure must be reduced by a process effectively lowering 
T < 4/3, (iii) a small initial m = 2 density perturba- 
tion must be present. The first condition may easily be 
met since collapsing primordial gas clouds are likely to 
carry substantial angular momentum. If an SMS forms, 
it is thus expected to rotate rapidly. The second con- 
dition is motivated by recent simulations of SMSs us- 
ing a microphysical EOS [6]. At sufficiently high central 
temperatures > 10 9 K encountered during the collapse, 
the effective T is lowered due to electron-positron pair- 
production. Even if the collapse is triggered purely grav- 
itationally, the temperature eventually increases, caus- 
ing an effective decrease in T. We find that lowering T 
by not more than 0.25% accelerates the collapse suffi- 
ciently to form a pair of two SMBHs. Finally, condition 
(iii) does not seem unlikely, since recent simulations of 
the direct collapse of primordial gas clouds [14 find that 
central supermassive objects that may lead to SMSs form 
m = 2 structures. In any case, small seed perturbations 
are likely to be present at the onset of collapse. For our 
SMS model, the m = 1 and m = 2 modes have the fastest 
and very similar growth rates [28, 29 . Therefore, even a 
tiny m = 2 bias in an initially randomly perturbed model 
may result in fast growth of the m = 2 mode and forma- 
tion of two fragments. A detailed analysis of this will be 
subject of future work. 

If m = 2 fragmentation and binary SMBH formation 
occurs in SMS collapse, the coalescence of the SMBH 
binary will result in a unique GW signal that can be de- 
tected at redshifts z > 10 with DECIGO and Big Bang 
Observer for SMSs in the mass range 10 4 — 10 6 M . If 
detected and identified, such a signal would confirm the 
existence of SMSs and could potentially inform us about 
their rotation and thermodynamics. Moreover, the in- 
teraction between matter and the binary SMBHs will 
likely result in an electromagnetic signature, the details 
of which will need to be established by future work. 
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FIG. 3: Convergence of the L2-norm of the Hamiltonian con- 
straint \\H\\ 2 of model M2G2 for the low rO, medium rl, and 
high resolution r2 (not rescaled). As the resolution is in- 
creased, the error decreases with a convergence rate between 
first and second-order as expected. The peaks occur when the 
black holes form, and when the black hole interior is not yet 
excluded from the calculation of the L2-norm. 
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Numerical Convergence 

We show numerical convergence of our results by con- 
sidering the Hamiltonian constraint H. H = is one 
of the four elliptic constraint equations which arise from 
the 3+1 ADM decomposition of the Einstein equations 
(e.g. [50 J. Any solution to the Einstein equation must 
satisfy these constraints. At the continuum level, if they 
are satisfied initially, they are automatically satisfied at 
all times. During a numerical evolution, however, nu- 
merical error can lead to constraint violations. As the 
resolution is increased, the constraints must converge to 
zero at the expected order of accuracy inherent to the 
numerical scheme. 

As detailed in [32], the overall accuracy of our sim- 
ulations is limited by the finite-volume hydrodynamics 
scheme, which is second-order accurate in regions where 
the fluid flow is smooth, and reduces to first order near 
shocks and discontinuities. Thus, we expect between 
first and second-order convergence. We note that the we 



TABLE II: Conversion factors for 1 solar mass Mq in code 
units c— G — M — K—W^o cgs units. 



Quantity 



code 



cgs 



Mass M 


M 


Length L 


M 


Time T 


M 


Density p 


M 



1.9891 x 1O 33 (M/M ) g 
1.4771 x 10 5 (M/Mq) cm 
4.9271 x 1O" 6 (M/M ) sec 
6.1716 x 1O 17 (M /M) 2 g cm- 



evolve the spacetime geometry with a fourth-order finite- 
difference scheme [32 , but the overall error is dominated 
by the fluid evolution. 

In Fig. [3j we show the L2-norm of the Hamiltonian 
constraint ||-H"||2 of model M2G2 (which forms a pair of 
black holes) for three resolutions labeled by rO, rl, and 
r2. The coarse resolution, rO, has 25% reduced reso- 
lution compared to baseline resolution rl, and r2 has 
25% increased resolution. The L2-norm is taken over the 
entire domain, excluding the interior of the black hole 
horizons. With increasing resolution, ||iJ||2 decreases 
with a convergence rate consistent with first and sec- 
ond order as expected. More specifically, we find that 
\\H\\ r 2 2 < 1.2511^11^ < 1.2511^115°. We note that the nu- 
merical evolution is initially non- convergent for the first 
t ^ 250M due to the non- constraint preserving applica- 
tion of an initial density perturbation. 



Conversion of Units 

For convenience, we provide the conversion factors 
from code units c = G = M = K = lto cgs units in 
Table \U\ As an example for a M = 1O 6 M star, a code 
unit time of T = \M corresponds to a physical time 
t ~ 4.93 sec, the initial stellar equatorial radius r e = 80M 
corresponds to c± 1.2 x 10 8 km ~ 0.8 AU, and the initial 
central density p c = 3.38 x 10" 6 M" 2 is ~ 2.1 gem" 3 . 
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